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S^aI^mary 


This is the final report on the theoretical studies of impact 
of composite plates completed by the principal investigator at Prince- 
ton University. Previous reports under this grant have presented analyses 
and computer codes for the calculation of stresses in composite plates due 
to central and edge Impact of hard objects. These studies were directed 
toward the problem of foreign object damage in Jet engine fan blades. The 
present report is directed toward three separate problems related to foreign 
object damage in composite plate like structures. These are; the effective- 
ness of constrained layer damping for leading edge impact protection; the 
development of multilayer mathematical models to calculate interlaminar 
stresses due to impact of composite plates; and a review of fluid modell- 
ing techniques for predicting impact stresses and forces due to bird im- 
pact. 

i “ Constrained Layer Damping of Impact Stresses, 
n a previous report an analytical-computational code was developed 
to predict the stresses due to the in plane edge impact of an anisotropic 
plate. This code included provision for an elastic protection strip to 
be placed between the impact force and the half plane of the plate. In 
the present report this code is modified to include a viscoelastic layer 
between the elastic protection strip and the composite plate. Similar 
techniques for damping plate vibrations have proved very successful. 

Numerical results show that a very thin elastomer damping layer may signi- 
ficantly reduce the normal impact stresses in the plate . The results are 
based on a modification of the plate-protection strip boundary condition. 
Since the code uses the fast Fourier transform, experimentally determined, 



frequency dependant material constants for the elastomer can be included. 

Part II - Multilayer Model for Impact of Composite Plates. 

In earlier studies by the principal investigator, the central impact 
of ccmposite plate was modelled using a plate theory that included linear 
bending and shear displacements and a single transverse displacement vari- 
able which effectively neglected wave propagation throu^ the thickness of 
the plate. In the present report hi^er order inertia variables are in- 
cluded. In addition the plate is broken down into c. set of identical 
orthotropic layers. Each layer may represent mauy plys or in specialized 
cases a single ply of the composite plate. Incorporation of these two 
features results in a model that can predict interlaminar shear and nor- 
mal stresses as well as wave propagation throu^ the thickness direction. 
Results for the line impact of a two layer plate show an interlaminar ten- 
sion developing under the load after impact. The computer code which 
solves finite difference equations for a periodic set of oscillators can 
handle any number of layers. 

Part III - Dynamics of Bird Impact. 

Prediction of impact stresses in composite fan blades not only depends on 
the structural modelling but on an accurate knowledge of the force be- 
tween the foreign object and the structure. In the final section of this 
report methods for predicting the forces generated during a bird impact 
with a solid wall are investigated. The physical properties of birds 
as related to impact are revievred. Simple Hertz impact calculations for 
bone and ccanposite materials show that the skeletal bones of birds will 
disintegrate under impact suggesting that a fluid model for impact, mip^t 
be useful for high speed impact greater than 50 m/s . A brief oiirvey of 
the literat\ire of rain drop Impact and computational fluid mechanics is 
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presented. A marker and cell hydrodynamic code is used to calculate 
pressures and force history between a flat plate and viscous incompres- 
sible fluid of cylindrical and spherical shape. Both normal and oblique 
impact are studied. The pressure and force histories show a fluctuation 
behavior suggesting either real or computational Instabilities in the 
code. Velocity distributions show the development of an eddy effect 
near the wall and a subsequent stationary zone of liquid near the plate. 
This leads to a fairly uniform pressure dist-'ibution across the contact 
area between the fluid and the wall. 
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PART I 


Edge Impact of a Plate With Constrained Layer Damping 


by 

F.C. Moon 






Part I - Edge Impact of a Plate With Constrained Layer Damping. 

To prevent failure of composite fan "blades under Impact forces, 
leading edge protection strips have been empltyed. In practice strips 
of stainless steel are wrapped around the leading edge of the blade. 

The effect of the strip is to "thwart the force of impact, "thereby spread- 
ing the normal impact stresses over a large area of the composite ma- 
terial underneath the strip. 

In a pre-vious repott we developed an analytical model for edge im- 
pact of a composite plate with edge protection [3]. In this model the com- 
posite plate was treated as a homogeneous anisotropic elastic material 
and the edge protection modelled as a beam which is connected to the edge. 
The results of "that study showed that the edge strip could decrease both 
the normal and edge wise impact stresses at the siirface but that signifi- 
cant shear stress could develop at the bond between the edge strip and 
the plate. The calculations were carried out for in-plane impact loads. 

The edge protection strip is only effective to the extent that it 
spreads the transient impact loading over a surface along the plate edge 
larger than the impact contact area. Also it disperses the pulse so that 
the Impact energy to the plate is spread out in time, hence decreasing 
the peak stresses. 

In contrast to the energy dispersion method of decreasing impact 
stresses the dissipation method uses damping material to convert the im- 
pact energy into heat instead of into kinetic and stored elastic energy 
in the plate. There are two approaches to the absorbtion of structural "vi- 
bration fa teams and plates. In one method a highly "viscoelastic material 
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is simply connected to one face of the beam or plate. Energy is converted 
to heat in the viscoelastic layer throu^ normal stresses. This method 
has been used to quiet the vibration of submarines. This method contrasts 
with the constrained layer method in which the energy is dissipated through 
shear stresses. This is accomplished by cementing a soft viscoelastic 
material to the plate and covering the damping layer with a stiff elastic 
material (see Figure 1 ). Thus the viscous layer ie constrained between 
two elastic plates, hence the name. In practice the viscoelastic sublayer 
is a thin hi^ damping elastomer while the constraining plate can be a 
thin plate of alumintm or steel. This method has been studied both ex- 
perimentally and analytically for vibratory motion by Yan [ 4 ] and Yan 
and Dcvell [ 5 ] • 

It is proposed to use such a method for the absorption of transient 
impact vibrations of composite fan blades by placing a thin elastomer 
material between the leading edge protection strip and the composite 
blade material 

In this section a model is proposed to examine the energy absorbing 
poteutial of such a concept for the edge impact of an anisotropic half space. 
The model is shown in Figure 2 . Between the anisotropic half space and 
the edge protection strip we asstune a thin layer of viscoelastic material 
which has a uniform normal strain in the direction £^3 aver- 
age shearing strain The inertia of the layer is neglected as well 

as bending moments. Thus the stresses t^^, and t^^ are transmitted 
from the beam to the half space unperturbed. However the campalibility 
conditions between the beam displacements U, W and the plate edge dis- 
placements ’o^, u^ are changed. 

If c is the uniform strain in the sublayer at position x^ and d 
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the thickness of this layer, then the relation between 


and W is given by 
(l-l) 


U3 - W . 633d 


Further if a is the change in angle of a line element normal to the 
damping layer (see Figure 2 ), then the relation between u^ and U is given by 



b ^ 

2 dx^ 


oc d 


( 1 - 2 ) 


However, the. angle a is not the total shearing strain. At the center 
of the damptog layer we take the average shearing strain to be 

’'13 - 

Finally to find 7^^ we use the viscoelastic constitutive rela- 
tions between the strains and ’^23'’ ^13 sublayer. We also neglect 

the strain in the sublayer. Tims for the Laplace transformed vari- 

ables ■vdiere s = ia> we have 

€33 = t33A(m) , ^3 . ya{c„) ( 1 . 4 ) 

where Y(oi), g(cd) are complex functions of the frequency oj, and t^^^ and 
tgg are neglected in (1-^). 

The stresses t^^ related to the svirface displacements 

of the half space u^, u^ through constitutive equations for the plate 
and the displacements must satisfy the wave equations in the plate. To 
find the displacements and stresses in the composite plate, we follow 
the same procedure as the no-stiip case except for the boundary condi- 
tions on the edge. In place of the zero stress conditions on the edge we 
8 
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relate the edge stresses t^^^ "to the motion of the beam strip, if 
one considers a small element of the beam-strip along the direction, 
the raomentTim balance equations in the x^, x^ directions become, (for a 
plate of unit thickness) 


Pb — s- = Eb + t, 


St" 


Sx, 
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(1-5) 


Pb 


St^ 


. El ^ + It, ..A 


h 

a sq «33 * 


In these equations U, V are the x^, x^ displacements of the beam 
element at the half thickness, and t^^, t^^ are the interface stresses. 

The compatibilicy condition between the beam and plate displacements 
W and u^ are given by 0^-1) and the condition between U, and u is 
given by (JC-^ . 

In the above equations b is the depth of the strip, E, I, I are 
lespecti/ely the Young's modiilus, moment of Inertia and rotary inertia. 

Also p^f(t)g(x^) is the edge loading applied to the outer protective 
strip surface. 

The equations for the plate remain as 'n the free edge case [ 3] 
and a solution is obtained by taking a Laplace transform on time and a 
Fourier transform on the space variable x^. With nondlmensionalization 
the solution in the plate has the form 


f-l 




r~ ”1 

“1 

i 

^ 1 

1 

■Vs , 


m * 
1 
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where p^, and are given by the following equa^xons (u 

indicates a Fourier and Laplace transformed \^riable). 


2, L> 


det = ^33^55^^ + tC23(-s'^-C^^kJ)+C^^(-s^-C^^H^)+k^(C^3+C55) ]p 


(1-7) 


+ (s^+C^^k^)(s tC^^k^) = 0 


where the C are the equivalent elastic constants for the anisotropic 

ii 

plate. We wxll choose the p’s with positive real parts to insure 
the decay in direction of the surface wave. Let the solutions be 

p = p^,p 2 ^ therefore, we have 


^ - C 


^(k^, s ) 





“ ’*^31^1 (l-8) 


P 



C2(ki,s) 



%2^2 


i[-s^-Ci^k^+C^^p|] (2) 


The equations of motion for the beam (l-5a,b), are next transformed 
using a Laplace transform on time and a Fourier Transform on the space 
variable x^. When the solutions for u^, u^, along with the constraint 
conditions for W, U are substituted into the bending and extensional 
equations of motion for the beam, and the constitutive equations (l^)^ 
used to eleminate the stresses t 23 j t^^ obtain two equations 
for the unknown constants C^, in (1-6) - These equations can be 
put into the form. 


r\ «ii 


pl1 






= P f(? 


^^2 ^2 - 
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To define the terras G,., we use the s mbol s for the 


Laplace transform variable and k for the Fourier transform variable 
where 

F(k,s) = / / F(x,t)e“®\t e"^^dx 

-00 o 

‘Then using the following symbols 
b^ = O.Ik^+I bk^s^-ipbs^} 

b/:^ = fEbk^+pbs^} , 

the matrix elements in (i-^ are given by 

G2(Pi,i|r3i) = b^(l+ 11^3^) - C^^(iki|r3^-p^) (I_10) 

+ b^d j^2Y(g ) (C33P^i|^3^-ikC^3)+ik(i^^^- ^j^[c_,ik-C33P^tY]) 

" G(s) ‘^55^^*3l''^l^J 

and where 

Hi = Oi(p2, (,32) 

”2 “ *32^ 

When the thickness of the viscous layer is sot equal to zero^ 
i.e. d ^ 0, then one obtains the solution for an anisotropic plate with 
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with a beam glued to the edge. When the beam thickness is set to zero, 
i.e. b s 0, then the free edge plate is obtained. These two cases were 
studied in Reference [ 3 1 • 

To obtain solutions in the time domain for a pulsed imput fg the 
expressions for u^, u^ must be inverted. Bais was accomplished using 
a double fast Fourier transform as described in a previous report [ 2 ] . 

In the example choosen a specific elastomer was chosen whose viscoelastic 
properties were known. The material choosen was an elastomer made by the 
Dupont Corp., LR3-6o4. This material was used by Yan in his disserta- 
tion [ 4 ] and the numerical values of Y(s) and G(s) were obtained 
from data in Yan's thesis. This data is shown in Figures 3,4. The shear 
modulus can be represented in the form 

G(s) = G(icu) = G'(cu) + i 

= G’(l+i G"/G') . 

The expression for G' (oi) can be represented by a cubic function of 
log(co/2jr). The ratio G"/G' is known as the loss tangent and for the 
particular temperatiire chosen can be approximated as a bilinear f'unction 
of log(o3/2jr). Thus for each frequency component in the Fourier inversion 
of the solution the corresponding value of G(s) was used. 

Results of calculations for a specific case are shown in Figures 5,6 
The plate material is a + 15 degree layup angle graphite /epoxy composite. 
The parabolic loading length is a = 2 cm and the steel beam strip thick- 
ness is b = 0.5 cm. The contact time in this (example is 35 dsec. Plot- 
ted in Figure 5 is the maximum normal stress at the plate interface 
versus the normalized thickness of the shear sublayer d/a. One can see 
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that ’While the str-ss rises for very small values of d/a, the noimal 
impact stress decreased dramatically for elastomer thicknesses less then 
20'^ of the impact half length. 

Figure 6 shows results for the edgewise stress t^^^. The maximum 
edgewise stress without the beam is about 4.8 from the study of Ref. 3 
and with the beam glued to the edge is 1.1 p^. Adding a damping sublayer 
appears to increase the stress for d/a <0.15 and the limit as d/a 0 
does not appear to result in the zero sublayer case. This is believed 
due to the fact that the shear modulus for the sublayer is many orders 
of magnitude below that of the composite or the beam. Thus the sublayer 
acts as a zero shear stress boundary condition . To check this we ran a case 
for d/a = 0 but the shear condition for the beam set to t =0 while 

13 

maintaining continuity of normal stress and displacement. This can be 
accomplished in the computer model by setting = 0 in equation (I- 10). 
The result of this calculation leads to a maximum t = 3.73 p I’or 

the same loading conditions and plate material as the cases above. This 
value appears to be the limiting value for d/a -» 0 and explains the 
apparent paradox. 
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Multi-Layer Model For Wave Propagation in Composite Plates Due to Impact 




F.C. Moon and B.S. Kim 
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Part II - Multi-Layer Model For Wave Propagation in Composite Plates 
Section 1. Introduction 

In our previous reports on wave propagation in composite plates [1], [ 2 ] 
the multiply plate was modeled by inclusion of linear bending and shear 
displacement across the thickness and a single transverse displacement 
variable for the midplane. This model is a modified Timoshenko plate 
using a procedure for obtaining approximate plate theories from the equa- 
tions of elasticity developed by Mindlin [6 ]. This simplified model as- 
sumes that the wavelengths of the impact forcing function are equal to or 
greater than the thickness of the plate. It is further limited in that it 
cannot treat wave propagation through the thickness of the plate and predict 
damage phenomena such as spalling. 

A number of researchers have presented models for a multi-layer com- 
pcisite plate. Many, however, have stopped short of the transient impact 
problem and have examined only the frequency-wavelength dispersion rela- 
tion for wave propagation in the p3ate [ 7 ]-[l2]. In this report we pre- 
sent another attempt to mathematically model the multilayer plate but will 
develop a method wherein propagation through the plate thickness can be handled and 
transier : impact stresses can be calculated using an inexpensive fast Fourier 
algorithm on the digital computer. 

The composite plate is represented by N layers; each layer may contain a 
number of plys (Fig.?). Each layer is treated as orthotropic with the symmetry 
axes of all the layers alligned. For alternating ply composites each layer 
should bontain two or more plys. The model can be extended to Include the 
case of the layer symmetry axes at angles to each other but will not be re- 
ported here. A key assumption is that all the layers are identical. While 
restricting the application, this assumption allows us to formulate the problem 
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using difference-differential equations. The tecshnique of periodic 
structures has been used in the study of electrical transmission lines 

p.3 ] and in the vibration of multistory buildings [l4] . A set of equa- 

■* 

tions of motion is developed for a typical layer. The relative motion 
of one layer to another is related by a phase shift. In this way the 
ninnber of layers can be increased without increasing the size of the 
matrices to be Inverted to satisfy the boundary conditions. 

The model incorporates the interlayer stresses as explicit variables. 
Throu^ these stresses we hope to extend the analysis to the study of im- 
pact of composite plates with viscoelastic damping layers and with cracks. 
Such studies are now underway. In the results presented in this report 
only the line impact has been treated. Ihis has simplified the calcula- 
tions and saved computer time in testing out the model. The technique 
however can be extended to the two dimensional or central impact problem. 
The next sections will describe the model in detail and discuss the numer- 
ical results. 


Section 2 . Formulation- 

Basic Theory of Linear Anisotropic Elasticity . Cauchy’s equations of motion 
in cartesian tensor form are 

^13, i " 

(ii-i) 

*ij ^ 

where body forces are neglected and the stress tensor is related to the 
infinitesimal strain tensor by 

+ C € 
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or in a condensed form this is often written as 


■ 


c ^ 

io i 


(n-3) 


elastic moduli 
materials 




or 


'ij 


has the following form for orthotropic 


"n 

^12 

=13 

0 

0 

0 

^12 

°22 

^23 

0 

0 

0 

“13 

=23 

Ss 

0 

0 

0 

0 

0 

0 

C44 

0 

0 

0 

0 

0 

0 

C55 

0 

0 

0 

0 

0 

0 

^66 


(II-U) 


Analysis of a Layer . For a layer shown in Fig. 7 we eB 5 >loy the approximate 
plate theory of Mindlin [6 ] and the displacement field u is expanded 
in terms of the Legendre polynomials as 


u(XiiXg,x^,t) = E u^^^(x^,X2,t)P^(ti) (II- 5 ) 

n=o 

where t) Is the local coordinate along thickness and nomeilized by b 
(b; a half of layer thickness). 

Instead of solving Eq.H-1 directly we obtain new approximate 
equations of motion by a variational process and integration over the 
thickness q. The result is 
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where 




M 


Idp (tj) 
J - 
-1 


dTi 




f « • 


By substituting the constitutive relation (n-2)with the displacement 
expansion (lI-5) into the above approximate equations of motion, we can 
find governing equations of motion in teiros of u^°^, u^°^, u|°\ ' 

The accuracy of the theory depends on how many terms of the dis- 
placement field we retain. Since the complexity in formulation increases 
rapidly with the number of terms included we keep terms only up to second 
order. Furthermore we will only examine harmonic waves propagating along 
the Xj direction so that we can drop terms and have [ } = 0. 

Next to get rid of the undesired coupling with hi^er modes we set 

( 2 ) ( 2 ) 

= lig =0. Then the resulting equations are 

2b(CiiU^^^+ + (■‘^21"‘‘^21^ = 2bpu^®^ 


^66 ^b ^,l'^^,ll^ = 2bpu^®^ 






( 1 ) 


( 1 ) 


^ ^ ^^12"^i, l"^ b^22^^2*^ ^ ) = 0 


(II-8) 
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Here we notice that the first, Iburth and. last equations are written 
in tenns of 'i \rtiere (n+i) is an odd integer eind represent the 
thickness stretching motion (or symmetric motion^ In the rest of the 
equations in iidiich (n+1) is an even integer, the u's represent flexual 
motion (or antisynnetrlc motion). Hence this process has decoupled the 
stretching motion from bending motion. 

To get rid of the 2nd order modes from Eq. (lI-8) we solve the 
last two eq\iations for and S'^d insert them into the remain- 

ing four eqiiatlons and drop the propagation of u^^^ along direc- 
tion which is equivalent to dropping 3^^^2l“*21^ equation. 

Then eq. (lI-8) can be reduced as follows: 


2b 


^66 ^b ^, 1^^,11 


+ (*22''*?2^ 



rii’^,! 



^^, 1 ^ 3 Cg 2^\2 *22 ^»1 


^*2l‘^*21^ 


(II- 



,(0) 




^* 22'*‘*22 ^ 




where 





Section 3. Wave Propagation 

Harmonic Waves . Let’s consider now a harmonic wave propagating in the 


direction. Namely one solutions for and . t are written as 


„(n) ^ 


i(kx,-<ot) 
t = T e ^ 


(II -10) 
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01. • ^isViL PAGE IS POOR 


In view of the Legendre polynooiical expansion the displacements on the 
both sides of a layer can be written as 


w- = ± } n = ± 3. 


+ ( 0 ) ( 1 ) 

= ;tj = ±1 


(Il-ll) 


and the displacement w and v can aI.so be given as 


w = W e 


i(kx^-ojt) 


V = V e 


i(kx^-cut) 


(11-12) 


If we substitute Equs- (11-10,11^2) into Eq. (lI-9) we find 


(-C^k^+5^)(w‘^-W“) + C^iK(v‘^-V") + b(T^-T][2) = 0 


-C^iK(w'^+w’) + i(-3C^^-Hi;2)(v’^-v’) + b(<^+T;^) 0 


3' -^ 22 


22 22 ' 


CggiK(W-"-w“) + (-CggK^-to2)(vV’') + b(T22-T"2) - 0 


(11-13) 


i(-C,,«2-3C..-f5®)(W^W-) - - ,iK(vV‘) +^iKb(T^,-T- ) 4 b(T^,+T- ) 


. 12 


3' ^11" 


o6 


3C, 


'22 


22 22 ' ''21 21 ' 


0 . 


^2 

where cS and k are defined by 


-2 ^2 2 
CD = pb CD 


K = bk 
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The ebove equation is the final form of governing equation written in terms 
of wave monber, frequency and anrplitudes of displacements on both sides of 

a layer. 

Plate Analysis . Remembering that above analysis is for any arbitrary 
layer in a plate, say the nth layer, equation (11-13) can be immediately 
written as a set of difference equations, [15l# 




=* 0 


= 0 

(11-15) 


= 0 


Vn-l^ ^'^n‘*^n-l^ " ° 


where we replaced bT^^ = t and bT^^ = c- 

tinuity conditions in displacements and stresses across the boundary be- 
tween layers are identically satisfied by these difference equations. 

Dlsnersion Relationship . For a plate made of R layers, in general, Eq. 
(11-15) gives 4 k equations written in terms cf U(N+1) variables 
T fO • . .,Tj^, Cjj)* Boundary conditions drop h variables among them 

so that Un unknowns can be determined by Uu homogeneous equations vhen 
the determinant of coefficient matrix vanishes, which provides the desired 
dispersion relationship between frequency cd, and wave length 2jr/k. 
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One layer Plate The dispersion relationship for a plate made of one 

Icyer (as shown in Fig. 8) can he found hy setting N « 1 in Eq. (15) with 

corresponding boundary conditions, i.e. c »t »cr, ct, eO and the 

o o 1 1 

resulting equations are now written in matrix form as follows: 



CigiK 

*1 

0 0 

\ 

r *) 

W-+W 
1 0 


0 

1 

O 

W 

|(-3C22-^^) 

0 0 

< 

Vt-V 
1 o 


0 

0 

0 



o 

H 


0 

^ 0 

0 

-CgglK . 




0 


(II -16) 


Then by setting the coefficient matrix zero we obtain the dispersion 
relations as 


- i(-CiiK2-3Cgg«;2)(. CggK^«=) = 0 


(II-17) 


Here we notice that the firs+ reodtionship corresponds to the state of 
deformation of and = -V^, which correspond to thickness 

extension of the plate or the symmetric mode, and the second describes 
the flexual defonnatlon or antisynanetric mode. The exact theory of plate 
gives an infinite number of dispersion relations but since we only 
kept inertia effects up to 1st order for b^th components of displacement, 
we have only the first four dispersion rela^J.on6hips . 

Dispersion relationships and corresponding phase velocities for an 
isotropic material with Poisson’s ratio l/4 (namely A = u) are given’ in 


Fig. 9a and b. Relationships for a 55^ graphite fiber-epo)(y matrix composite 
with layup angle 0* and 45" are shown in Pig. 10a, b and 11a, b. Note here 
that midifled frequency o5 and phase velocity c are normalized by 

Cgg. From these figures we clearly see that /k (or c^) approachs the 

limit which is the dilatation wave speed in case of an Isotropic 

plate and quasi-dilatation for anisotropic plate when the wave number « 
becomes large (or the wave length becomes small compered with the layer 
thickness b). Also notice that o3j^ /k (or c^), which for kb « 1 is 
a bonding wave approaches a shear wave for kb » 1 . Pig. 11a and b show 
the effect of anisbtropy on the dispersion relationships and the wave speeds 
of the dilatation (and quasidilatition) waves when the layup angle changes 
from O' to 90" (see Table II). 

* lu this case we obtain 8 equations by putting n = 0 and 

1 in (II-15). The boundary conditions require t = 0 = t = 0 = 0 fsee 

0022 

Figure 8). The 8 equations are written for 8 unknowns (W ,V ,W, , V, , 0 , ,w .V ') 

' o' o' 1' 1' 1' 1' 2 * 2 ' 

and again by following the same procedure as in one-layer case we find the 
dispersion relations as 


i { (-C3^K^-3Cgg-to2)+ 


" i £^66 
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r-c^ 

^ ^12 


p p 

K +(W -C 


'jj ) (“( -SCgg-to" )+ (Ji=='-CggK‘=' ) ) ] 


-2> 


^T.2 
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) - (» -C66K®)(-3C22«^)f(0>=-Cj^K2). 


+ . 0 


Dispersion relationships for the isotropic plate and an-’ ropic plate 
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with layup angle 0" and hT a.re plotted in Fig.l3a,h,l4a,b, and 15a,h. 


K-Layer Plate . In general, we can otitain a 2(N+1) order polynomial of 
U3 by expanding a (4n) x (4 n) determinant and find 2(lH-l) dispersion 
relations. But unfortunately this process involves consideraoly compli-> 
cated algebra and it may be necessary to develop a computer technique to 
find roots of an equation in determinant form (not in polyncanial form) . 

A difference equation approach can be used to solve the N set of 
four simultaneous first order difference equation given by Eq. (11-15). 

This proceedure is neat and can be generalized for any ntimber of layers 
but the last step of this approach, where a long polynomial is to be solved 
again, is not simpler than previous direct method. 


Section k-- Difference Equation Approach for Inrpact Problem 
Solutions of Differ ence Equations . Since the simultanious difference equa- 
tions given by Eq. (II-I5) are linear and all the coefficients are constants, 
we can try the following form for their solution [15} 


T = A e‘ 
n 


2ipn 


"n = Se 


2ipn 


(11-19) 


W 


n 


C e 


21Sn 




2i3n 


where 0 is complex, in general. By substitution of these solutions 
into Eq. ( 11-15) we have 

r 

. . p -p. 

-C^Ksin0, 0 1 sin0 


2 ^2 

(-C^K -KD )cOS0, 


2 - 2 . 


(“Cgg ; +<^ )cos0, i sin3 


3 ( )sin 0, -CggiKCosP, - Ksin0, cos0 

^^22 

- C^iKCose, |(-3C22+w^)sin0, cosS, 0 


C 


D 




^ '1 

0 


0 




ej 


0 

V. j 


(11-20) 


Again if we set the value of the determinant zero we find the follow- 
ing equations for B, ui and k • 
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( 11 - 21 ) 


k 2 2 k 

a^cos p + 3 cos 3 + a^sin 3 = 


/- 2-2x/ 2-2v 

&i = -HU )(-CggK •Kd'") 


^2 

“3 =%-3<:22«^) 


= (-C^A i(-0 K^4<52)(-3a„„-M2)) 




^ flf ^ 2 ,-2s. ^12^66 2, 

°3 - *'^^ 66 ’^ ^ ^ 


3C, 
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“1, = (=66*(-%6"^“'>3cr 
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Then we can find k values of 3(say ± 3^ (±3) and ±3g(± a)) with given 
values of ai and k. Accordingly solutions given by (11-19) can be writ- 
ten as 

+ A + A + A,e-2i“ 


"n = \ 


2i3n . „ _-2i6n . „ _21ocn . „ .-2ic£n 


e + B^e 


+ B^e + Bj^e 


(11-22) 


«n ' + CjB^K* + Ci,b-2‘“" 
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Next Tflhen we substitute these solutions to our original difference equations we 

can find the relationships between and • The boundary 

conditions on the top and bottom faces of the plate reqijire calculation 

of T and a , and they are given bj 
n n' 






„ -21^n . ^ 2iocn ^ „ ^-Siotn 

Bge ^ + B^e + Bj^e 


with 

-Ko^ ) cos^p-Kx^iJ^s inS 
= (ctgC^ m k; cosPsinB 


(11-23) 


(11-24) 


vhere the unknovm constants B^^'s have to be determined from boundary 
conditions . 

Dispersion Relationship and Ijnpact Problems . The dispersion relation for 
a composite plate consisting of N layers can be found immediately by 
setting ■’’o “ ~ ^ which leads us to 


1 

0 

1 

0 

0 

X(P) 

0 

X(a) 

cos2PN 

1 sin2PN 

cos20W 

i sin2ow 

X(P)sin23N, 

X(P)cos28N, 

i X(a)sin2Q!N, 

x(a)cos2cilN 


= 2X(a)x(e)(l-cos20W coB2BN)-[5^(a)+X^(P)}sin2Pn sin20W (II-25) 

= 0 . 
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where P and a are obtained by solving ( 11-21). 

This difference equation method can be applied to the impact problem 

without any further difficulties by using integral transforms (Fourier 

transform in and Laplace transform in time) instead of harmonic 

wave analysis. The resulting eqmtions are the same as (11-23) and B^'s 

can be determined from a . t * a„, which are now the integral trans- 

o o' W N 

forms of impact functions. For the present report we have chosen a line 
impact along the x^-axis, i.e., 

t„ = -P (l-(“) }sin -2^. ; on + side of layer, (II-26) 

22 o a 

for -a < ^ 

for the only nonvaDlshing Impact function. Therefore the resulting 

A 

boundary conditions are ^ “ ^22 

integral transform of t^g* Once the B^*s are determined, the dis- 
placement fields and stress fields can be computed by inversion of the 
integral transform. For the present problem the inversion cannot be 
accomplished analytically because of the complexity of transformed func- 
tion, but since the impact function given by Eq. (11-26) has finite rise 
time, d\iration and extent both in time and space, inversion can be carried 
out numerically by use of Fast Fourier Transform techniques. 

Section 5» Numerical Results 

Numerical inversion of the solution for the stresses in a multilayer 
plate was carried out for a two layer model of a composite plate. Each 
layer may contain many plys, but for the unidirectional fiber layup case 
each layer may represent a single ply. The calculations were carried us- 
ing equivalent anisotropic elastic constants for a 55% graphite fiber/epoxy 


matrix composite plate. A two layer model allows direct calctilation of 
midplane interlaminar shear and normal stresses. 

The propagation of a wave aftej, impact on a plate consisting of two 
steel ].ayers is shown in Fig. l6 a ~ f where we can see two distinct 
wave speeds; c^(= 5-33 mi/\isec) for w's and cr^ and Cg(- 2.67 mm/|isec) for 
v's and while the dilatation wave in steel has a speed of = v(X+2n)/p st 
5.61 mm/^sec and the shear wave c^ s/p/p = 3.25 mm/psec. This Indicates 
that the initial signals are propagating via the acoustic "branch of the 
symmetric mode with an almost constant group velocity c. = 1.63 c = 5*31 

«L S 

nm/psec and the major signals are carried "by the bending mode (the acoustic 
branch of the antisymmetric mode) whose group velocity is lower than c 

s 

(as shown in Figure 9). Similar phenomena is also observed in case of 
an anisotropic composite. 

Figure 17a, b shew the interlaminar shear stress versus time for the 
45“ fiber layup cese (load perpendicular to the fibers) and the change of 
interlaminar shear stress along the plate at various times after impact. 
Figures l8a, b present simular nmerical data for the Interlaminar normal 
stress. In Figure l8b one can see that directly under the load the normal 
stress is initially compressive but subsequently becomes tensile. This is 
due to reflection from the back surface which in the two layer plate model 
is an oscillation in the thickness direction. Such tensile stresses may 
account for spalling damage and ply separation. 

Finally in figure I9 data are presented for the case of the load in 
the direction of the fibers. Here for the case of Interl am inar shear t^ 

Y 

Propagations of w^, , v^ and v^ are almost exactly same as 

those of Wg, and v with different signs in case of u and 

xj-o d 

and they are not shown here. 
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one can see a distinct wave spreading along the plate in time. From 


the figure we find c = I.3 mm/^sec which is slightly lower than \/Cgg/p 
as in the isotropic case mentioned before. Investigation of wave propaga- 
tion throu^ the thickness direction requires an Increase in the number of 
layers and is underway at the writing of this report. 
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Part III - Dynandcs of Bird Impact with Aircraft Engines 
Section 1. Introduction 

It is well known that hird impact on the fan blades of Jet aircraft 
poses a serious 'threat to airline safety. Ihis problem has been studied 
extensi'vely using both dummy birds and actual carcasses, both in Great 
Britain [16] and in this country Dl-TU^It has been clearly demonstrated 
in these tests that gross damage to composite blades can occur on impact 
producing broken parts of blades which themselves can initiate sequential 
fracture of the rest of the blade set. Films of single blade en- 

counter with bird carcasses or simulated bird materita suggest that the 
bird may be modelled as a transient fluid mechanics problem. However, 

■until recently very little analytical or computer modelling of bird 
impact was available in the tediniceil literature. A large literature on 
computer modelling of transient fluid mechanics prb.lems exists with ap- 
plication to rain impact and erosion [19-22] but little if any had 
been applied to the bird impact problem. 

The objective of the Princeton program in this area was to search 
the fl'uid mechanics literature for solutions and computational techniques 
•that could be used to predict the forces and pressures on the blade struc- 
ture during bird Impact. We had also hoped to use such forces to simul'tane- 
ously predict both fluid (bird) and blade motion (and hense stresses) dur- 
ing impact. These goals were only partially met as will be discussed be- 
low. But the principal problem lies in the reliability of the forces and 
pressiires obtained from the ccmputer simulation programs. 

Before a proper model can be chosen, one must examine some of the 
physical properties of ccmmon birds. A complete description would include 
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the effects of bones and feathers and noncalcified tissue such as muscle, 
tendon, and fat a? well as vital organs. A mechanics description of such 
an object wc^uld include such descriptions as inhomgeneous, viscoelastic 
and nonlinear. A complete solution of the impact of such a material is 
not possible at this time. Using contempcrary techniques, one can hope 
to obtain a fluid model which is homogeneous, viscous, and perhaps com- 
pressible. 
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Section 2 . Ebyslcal Propertlee of Birds 

Before examining potentisG. mathematical models for bird impact 

studies we rewiew some of the physical properties of common birds* A 
summary of wel^t and geometric properties of birds was given by (^ifflths 
[23]* Elastic and inelastic properties of bones and muscles of many ani- 
mals and birds has been compiled by Yamada [2U] . Ultrasonic wave pro- 
perties of fatty tissue and muscle material are fo\ir*l in a review by Fry 
and Dunn [25]* For further details the reader is directed to the growing 
literature in biomechanics, in particular the collection of reviews 
edited by Fung [26] . 

A stnnnary of the information found in these references is presented 
in Table 1 . It should be cautioned that the numerical values given are 
in general r»ugh values sind in some cases may not be representative of a 
class of birds because of the small number of specimens sometimes tested. 

In summary xhe wei^ts of birds range from 9 (20 lb • ) for a swan 

to 1/4 k0n (0.55 lb.) for a sparrow hawk. While the density of mammalian 
fat and muscle is close to that of water the overall density calculated 
by Griffiths [23] was found to be less than that of water. He attributed 
this to air sacs which he has estimated range from 10-20^ of the volume 
of pigeons and ducks. Hie skeletal structure of flying birds comprises 
less then 10^ of the wei^t [23]. More extensive data for chickens may 
be found in [16], since these are readily available in the commer- 

cial food industry. However, data based on chickens which are ground 
birds may be misleading if extrapolated to flying birds, which are often 
involved in foreign object damage to aircraft. 

Hie ultimate compressive strength of the femur bone of domestic 
fowls and birds is about 686O N/cm^ (9,950 psl) in the longitudinal direction 
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and atout 45^ less in the transverse direction (Tables 26, 30 Yamada [24]). 

The elastic modulus in compression for an ostrich femur is 0.5 MR/cm^ (0.715 10 psi) 
[24]^ while the value in tension is 1.3^ MN/cm.^(2.0 10 psi) in the wet 
condition. Values for other birds were not available. 

The strength of bones under Impact may be enhanced by the protec- 
tion of skii. . Currey [ 27] found that 37^ more energy was required to 
break rabbit bones protected with skin under impact than those without 
protection. 

While bone may be treated as an elastic material, muscle 
and tendon are hi^ly nonlinear materials. The ultimate tensile strength 
of tendon for domestic ducks is around 6370 N/cm^(9200 psi) with 6.7^ 
elongation (Table 73 Yamada [24]). 

The ultrasonic wave speed in mammalian fat and muscle is around 
1500 m/s whidi is near that of water [ 25] . However the decay of ultra- 
sonic waves in fat and muscle is much greate* . At 1 MHz the characteristic de- 
cay distance is 7 and 20 cm for muscle and fat respectively compared to 
4000 cm for water. Thus the water hammer model , employed in rain impact 
problems, which has a shock wave generated in the watf-r on impact, may 
not be approp riate for bird impact because the large damping would smooth 
out or imnede the attempt of the waves to form shocks . Further the pre- 
sence of bone would further disperse any shocks developed by scattering 
the waves. 

The viscous nature of soft tissue is also much greater thari water and 
is estimated to be as hi^ as 150 poise compared to 10 ^ poise for that of 
water or 15 p>oise in the case of glycerine. However a more realistic model 
would certainly include viscoelastic effects which have been measured for 

certain biological materials [ 26] pnt are not renor^ -d here since only fluid 
•models for birds will be discussed. 
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Section 3 . Bone Impact Model 

If the bird is to be treated as a liquid it must be shown that the 
impact of skeletal structure of the bird with the fan blade will gener- 
ate stresses greater than the strength of the bone material. To model 
this we consider the bone as an elastic cylinder of radius Rg under im- 
pact with the fan blade material of radius of curvature R^. Treating 
the bone as elastic will obtain an upper bound on the stresses that would 
have to be sustained by the bone in order to remain intact. 

The solution for the impact of two cylinders, as shown in Figure 20 
may b^ fo\md in the monograph on impact by Goldsmith [ 28] . Kie impact 
theory presented in [ 28] is based on that of Hertz \rtiich starts with 
bhe contact force between two elastic solids 

where ot is the relative approach of the two bodies and kg is a con- 
stant depending on the elastic constants of the compostt.' end bone, and 
the geometry Rg* results of this theory is the time of 

contact T 



where V is the normal velocity and M is the mass of the bone cylinder. 

The results are shown in Figure 21a,b.0ne can see that the contact times 

are less than 10"^sec compared with the time of flight of the bird mass 

-3 

past the blade of around 10 sec. 
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The maxlmnn ccmpressive stress can also be calculated with the Hertz 
theory. !Hie details are contained in [28] and are not shown here. The 
results are shown in Figure 21b for two cylinders with their axes at 90* 

to each other. Here we used = 0.64 cm for the radius of the hone, and 

P 

20 gn for the mass of impacting cylinder and 0.4l MN/cm (.6 x 10 psi) 
for compressive elastic moduli of the bone [24]. The blade material was 
assumed to be graphite/epoxy. As one can see, the induced 
stress is order of 10^ Newton/m , (lO-"^ psi ) when the impact velocity is 
somewhez’e around 100 m/sec. This stress is much hi^er than the ultimate 
strength of the bone in a transverse ccmpressive load (5*10^ Newton/cm^) 
which implies that the failure of the bone is immediate. 

This rough calculation supports the idea that at speeds greater than 
50 m/ sec Uie bird may be modelled as a fluid since in any encoimter of bone 
with the blade the strength of the bone will be greatly exceeded. 

Section 4. Liquid Impact Models 

Bie Impact of a liquid object with a solid target has been the sub- 
ject of study in problems of rain erosion [29] -[34] and micrometeorite 
impact in the hi^ speed limit where the impacting object can be 
treated as a liquid. While water is usually treated as a nonviscous 
incompressible ixuid at low speeds, during the hi^ speed impact of 
rain drops the compressibility of the fluid becomes important and a 
shock wave propagates into the fluid upon impact with the solid in a man- 
ner similar to vaterhammcr in a pipe, Figui’e 22. If the impacting 
fluid is moving with velocity with respect to a rigid target, a one 
dimensional analysis of this problem predicts a pressure p given by 

p = P V V 
^ o o s 
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where V is the velocity of propagation of the shock wave in the fluid, 
s 

and Pq is the density in the \incompressed part of the fluid. If the 
target is elastic the normal pressure on the solid is [34] 


where 


P = 


p V V 
*^o o s 

P V 
e e 


D is the density of the elastic medium and V is the 
■the speed of sound in "the elastic solid* 


Biis analysis neglects motion of the fliiid lateral to the incoming 
velocity. In fact if the speed is low enou^ the fluid will flow tangen' 
tial to the surface rather than ccmpress noimal to the surface. In this 
hydrodyiamic limit the maxirium pressure is proportional to 



If enough fluid is present and some quasi steady flow is established 
near the center of impact a stagnation flow solution found in many bocks 
in fluid mechanics can be used [35], [36]. In this solution the tangentid.1 
velocity along the solid increases linearly with distance from the center 
of impact. 

Taylor [ 36 ] has shown that for the steady flow of a two dimensional 
jet of IncompresBible invlsid fluid, impinging on a flat plate, the maxi- 
mum pressure is pV^/2 and occurs at the center of impact (Figure 23 )• 

The pressure falls off by about 75^ at a distance equal to the width of 
the jet. Taylor has also presented data for oblique flow of a jet over 

a plate [ 36 ] . 
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A number of computational models have been proposed to solve the 
equations of fluid mechanics for transient problems. A review of these 
techniques is given by Roache [ 37] • Amsden, Harlow and coworkers have 
developed an extensive computational scheme to solve transient incompres- 
sible viscous flow problems [19]-[2l] ^ [38] as well as for compressible flow 
prob3.ems. In one published example they have treated the transient splash 
of a liquid drop into a pool of water as well as the rigid plate impact. 

In these examples they have neglected viscosity. Their results show a 
radial velocity increasing linearly with radius similar to steady two 
and three dimensional stagnation flow. 

Recently Huang, Hammit and Yang have presented a numerical scheme for 

a nonviscous ccmpressible fluid and have published the results for the 

impact of a liquid drop onto a plate [22]* The solution predicts wave 

propagation into the liquid but no propagation into the solid target and 

no shocks or surfaces of velocity jumps are included. The results are 

quite extensive. However in the published discussion following the paper, 

Ilejanann disputes the findings, claiming that shocks should be foimed and 

that waterhannner pressures pV V should be reached. In [22] the cal- 

o s 

culated pressures in the fluid and on the plate ere far below the 

theoretical waterhammer pressure. Also this program does not include the 
effects of viscosity which mi^t be important to bird modelling. 

Experimental studies of liquid impact pressures have been made including Brunton 
[31] and Smith and Kinslow [39]. Experimental bird simulation experi- 
ments have been preformed by Allcock and Collin [16] in Great Britain in 
which they meastire the force history, ^ftiey find that the maxiTmim force is 
proportional to the kinetic energy of the bird or square of the initial 


velocity. Similar results have been reported by Hopkins in the United 
States [ l<.o] . 
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The depeudance of impact force on the initial kinetic energy of the 
bird can he explained using a rou^ momentum analysis. Thus if is 

the average Impact force, at the time of impact, and all the momentum 
under normal impact is turned 90® to the initial velocity vector then, 

F^At = MV^ . 

If we choose the time of fli^t d/v^ for the impact time At , where D 
is the diameter of the spherical bird say, then 

M V® 
o 


This model can be refined a little by assuming that the momentTmi is 
changed during impact at a rate proportional to the rate at which the 
bird volume crosses an immaginary plane surface. Thus if the bird is re 
presented by a ellipsoid with a surface given by (see Figure 24) 



1 


then if z is the distance along the symmetry axis of the ellipsoid 
from the impacting tip of the moving liquid the impact force is given 
ty 

F = Vp nr^(z)^ = Pjt r^(z)V^ 


where r is given by the previous equation. 

For a sphere the maximiBi force is given by 


max 


2 D 
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o p 

At c 100 m/sec, D * 8 cm end p = lO’^ kgm/nr (density of water), 

^max ~ Newtons (- 10 Ibft) The average pressure of suck a force 

over an area equal to jtD A would be = pV^ = 10' N/m (1450 psi). 

!Qie waterhainmer pressure assuming a shock wave formed in the water 

Q p 

(Vg = 1500 m/s), would be p^ = pV^V^ = 1.5 10 N/m (21,800 psi). 

From the experiments of Alcoch and Collin [l£] and Hopkins [4o], 
the dependence of Impact force on would imply that the average pres- 
sures were also so dependent and that the incompressible model would be 
appropriate for birds . The compressible model with waterhammer pressures 
would lead to a linear dependence cf force on velocity. 

However, since the discrepancy between the incompressible and com- 
pressible pressures are so great, further study bn the effects of com- 
pressibility would be worthwhile. If the results of Huang, Hammitt and 
Yang [22] are proved rif^t - namely that compressibility does not re- 
quire shocks in the impacting liquid - then the experimental resiilts on 
bird Impact [l6], [40] might be compatible with a compressible model. 

Section 5. Results of Hydrodjh^amlc Computer Model 


The equations of fluid mechanics were solved using a finite differ- 
ence technique for both plane and axisymmetric motions. The differen- 
tial equations for incamprcss^.ble viscous flow are given below (see e.g. 
[35] and Figure 25) 


^ 3^ Sr^^ u^ ^ 5uv 

^ a "br dz 
r 



+ V 


bibu bVs 


dv . 1 ^r uv . dv 
St * 

r 



} 


u and V are velocities in the r, z directions respectively, 
g is the body force per unit mass and q!i the pressure/density ratio, 
a = 0 in plane coordinates; a = 1 for axi-symmetrlc flow or cylindrical 
coordinates . 

Incompressible flow requires 


D 


1 5r°^ u Sv 

a ^'r 5z 

r 
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The finite difference algorithm used to sdmulate bird impact was 
Amsden and Harlow’s simplified Marker and Cell program (SMAC) [20] with 
modification of the plotting routines, and boundary and initial conditions. 
The program was modified to accomodate the IBM360-91 computer and associated 
output devices including printer, plotting and microfilm hardware. 

Initially 50 x 50 or 50 x 30 cells were set up. Each cell contains 
nine marker particles. When the number of marker particles per cell is 
less than 9 the cell is designated a surface cell. The fluid "bird" 

occupied up to 300 cells. 

Both diffusive eind convective sources of numerical instabilities in 
finite difference methods require that the time and distance intervals 
at, As satisfy certain ineqvialities as necessary conditions for stability 
[21], i.e.. 


vAt 

(Az)2 


1 

2 


V At 
o 

Ax 


< 1 
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uhere v la the )cinematic viscosity and the initial velocity. A 

discussion of the stability of the MAC method has been given by Doly and Procht [4l] 

In calculations the above inequalities were satisfied. However 

a long time behavior of seme of the numerical results did not always 

eudilbit continuity in either position, velocity or pressure from time 

cycle to cycle raising questions about the reliability of the method. 

The SMAC method allows the use of either free slip or no slip boundary 
conditions, both of which were tried. While the bird is highly viscous, 
the skin and feathers mi^t provide an effective free slip boundary 
condition . 

Viscosity was kept as a parameter in these studies which was ignored 
in the splashing drop paper of Harlow and Shannon [191 and Huang et al. 

[22]. If D represents the diameter of the fliiid cylinder or sphere, 
the lniti€il velocity, and v the kinematic viscosity, then the 
Reynolds number 
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used in the computer simulations ranged from 10 < R < 10 . As a 
example we used the data, ■ 100 m/s, v » 9*5 cm /sec 

h 

viscosity of glycerine, density of water), D = l6 cm, R = 1.7 10 

2-5 

and At a 5 lO"® s, Az = 1 cm. For this case vAt/Az * 4.8 10 

and V At/Az = 5 lO"^ which are well below the stability criteria* 

o 

Two geometric configurations of fluid and target were studied. In 
the first, normal impact was studied for a fluid cylinder or fluid sphere. 

This geometry requires a solution for only half the fluid slug because of 
the Inherent syiaaetry in the problem. Fig. 26 shows the time sequence of fluid 
and surface cells of a half sphere under normal Impact with a rigid wall. In 
the second configuration a rigid rectangular target was set up and an cylindrical 


fluid slue oould Uspsct «ie target at either normal or ohUque incidence, 
in each case the fluid has an. initial condition of unlfom velocity with 


gravity ignored. 

Figures 27, 28 shows a time sequence of marker particles for normal 
ijnpact of fluid ctJlndar. -Jlie marker particles are shown in Figure 2? 
and the velocity vectors shown in Figure 28. Figure 29 shows oblique 
impact of a fluid cylinder on the edge of a rectangular object. 

A time sequence of marker particles and velocity vectors is shown in Elg. 
ures 30 , 3 Lfor normal impact of a spherical fluid slug with a rigid plate. 

For early time after impact the radial velocity at the plate shows 
a linear increase with radius which is charactersitic of stagnation point 
flow [ 35] . However as the impact proceeds there appears to develop an 
eddy current flow near the plate creating a dead zone of flT^d. This 
can be seen in the velocity plot in Figure 32 , and the radial velocity 
plot verses radius in Figure 33 • Thus if the eddy flow is physical and 
not due to numerical instability, the normal impact velocity in the fli^ 
actually reverses . A plot of normal velocity flow versus distance along 


the y axis for various times is shown in Figure 34. The velocity 
starts out unifoma and then the normal velocity of the fluid near the 
pinte approaches zero for small times and fin^y reverses flow for latet 
times indicating an eddy flow. This in effect produces a rounded station- 
ary fluid obstacle which deflects the remaining fluid away from the cen- 
tral plate impact point. This stationary central zone then tends to 
create a pressure that is fairly uniform with radius. 

The stress in the fluid is given by 



where (Uj^3 are the cartesian components of the velocity vector, fi is 
the viscosity coefficient, and p is the hydrostatic press\ire. The 
pressure p in the SMAC finite difference scheme is found from an itera- 
tive procedure and is a direct output of the program along with the 
velocity vectors in each cell. A plot of pressure to density ratio 
versus radius is shown in Figure ^ for a Reynolds number equvalent 
to a l8 cm diameter fluid sphere moving at 100 m/s with the viscosity of 
glycerine at various times during impact up to about 1.0 ms compared to 
a time of fli^t of 1.8 ms. The pressure versus radius exhibits fairly 
smooth behavior for a given time, and somewhat constant pressure versus 
radius for time between 200 |is and 800 ps which was suggested by the 
eddy flow phenomena. However the center pressure versus t;Lme does not 
show a smooth behavior, at first increasing then decreasing and finally 
increasing again implying a hi^ total force at the end of impact than 
at the beginning. 

Since we had intended to used the total force to calculate the rigid 

body motion of the target (fan blade) we attempted to check the computer 
calculated pressures and resulting force using a different technique such 
as integrating Bemoullis equation for the pressure. 


This equation Involves calculating accelerations 8v/8t which must be 
found from two sequential time solutions for v. The accelerations 
calculated in this manner however were not reliable and did not lead 
to a check of the pressure distributions. 

Another attempt involved adding up the total momentum H!v(l,J) 
overall the cells (Figure 36 ) , The total normal force is then 


> 


F = p g^\l,J)-v<^^(I.J)l _ 

n tg-t^ 

This too produced an erratic force behaTior and did not provide a satis- 
factory way to check the calculated pressures. 

However if a smooth curve is fitted to the momentum vs. time data^ 
(Figure 36 ) and the force calculated from this function, a continuous im- 
pact force history is obtained. In addition this force history ccanpares 
reasonably well with the integrated pressure profiles found from 
the numerical calc\ilation (Figures 37, 38). The force at first peaks 
and then attains a constant value for times up to about 20^ of the transit 
time of the fluid cylinder. Thus while the pressure-time data from the 
finite difference code is erratic from cycle to cycle, it appears to be 
at least consistent with the velocity or momentum data when averaged over 
a mmiber of cycles. 

The effect of slip or no slip boundary conditions on the pressure 
distribution on the plate is shewn in Figure 39, for a fliiid cyclinder under 
normal impact. For early times the pressures are about equal but beyond 
200 qsec the free slip impact results in higher pressures. 

Another observation for the full cylinder case is the development 
of unsymmetrical radial flow along the plate for normal Impact Figure 40. 
(such symmetry is of course guai*anteed for the half cylinder or sphere case). 
While such instabilities may develop in an actual flow, in the numerical 
solution this unsyrametrical flow probably indicates a numerical instability 
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In the finite difference code. 

In smunary velocity plots of viscous incompressible fluid Utrpact with 
rigid obstacles using finite difference codes would appear to offer a 
way to calculate the forces due to bird impact on fan blades. However 
lack of any exact analytical results to check the calculated pressures 
and forces raises doubts about the efficacy of using this approach to 
predict deformation of fan blades. The experience of the rain ijnpact 
problem,ln which there is great controversy over the actual pressures 
produced during impact, suggests that finite difference codes may not pro- 
riie a defiiiitL^<B answer for the bird Impact problem either without 
farther analytical, experimental or other computational check such as 

a finite element analysis. 
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Conclusions and Ke commendations 


1. The analytic modellin/r of constrained layer damping as a mechanism 
for decreasing stresses in composite plates due to impact shows promise 
of significant reduction of stresses for edge impact forces. It is 
recommended that constrained layer damping he studied for central impact 
of composite plates. To test these results, it is suggested that 

a limited experimental program be initiated on the concept of shear 
layer damping of impact stresses. 

2. The multi-layer generalization of a Mindlin plate appears to be a 
straif^t forward method of modelling a multiply composite plate for the 
study of impact response. The combined use of finite difference techniques 
in the thickness direction and the fast Fourier transform in the plane of 
the plate results in a fairly efficient method of studying interlaminar 
stresses and wave propagation through the plate. This technique might 

be modified to investigate the effect of interlaminar cracks or flaws on 
the impact stresses in the plate. 

3 . Films and calculations of stresses in bird bones due to impact seem 
to suggest a fluid model for the study of forces duo to bird imnact of 
aircraft structures. However analytical solutions for transient impact 
of a slug of fluid are not known. As shown in this report finite dif- 
ference comnuter codes can be used to obtain velocity, pressure and 
force histories. These "computer experiments" show the development of 
instabilities and eddy flow, in the fluid during impact. Wliether such 
motions are real or due to computational instability cannot be decided 
without comparison with either experimental results or other numerical 

h9 


schemes such as the finite element method. 

A search of the computational fluid mechanics literature reveals a 
her of potentially useful computer codes for the study of bird impact 
forces. These codes, if they proved accurate, caold save considerable 
sums in experimental testing. However several questions concerning 
numerical stability and accuracy of the impact pressures and forces must 
be carefully examined before they are embraced. While a bird is made up 
of highly viscous materials, the effect of compressibUty and of shock 
propagation into the fluid bird needs to be examined. 
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TABLE I 


Riyslcal Properties of Birds 


Name 


Weight (kgm) 


'Reference 


Caramon Gull 


0.45 


[23] 


Wood Pigeon 


0.45 


[23] 


Mallard Duck 


1.1 


[23] 


Canada Goose 


3.8 - 6.4 


[23} 


Whooper Swan 


9 


[23] 


percent wei^t for e dcken 

Body 

Wings 

Legs 

Head + Neck 

Ref. 


67.0 

6.4 

8.6 

8.0 

[16] 




Specific Properties 



Name 

Density 

Speed of Sound 

Rate of Decay 

Viscosity 

Ref. 


gjn/cm^ 

m/s 

(cm“^) 

•poise 






(lO’^S/m®) 


Mammalian Tissue 

1.07 

1570 

0.13 (@ 1 MHZ) 

150 p 

[25] 

Human Skull bone 

1.7 

3400 

1.7(@ 1.2 MHZ) 
7.8(@ 3.5 MHZ) 

chicken blood 


Chicken Body- 

1.05 



3-5 10-2p 

[16] 

Glycerine (20®c) 

1.26 



15P 

[25] 

Water (20® C) 

1.0 

1500 

25 10 ^ 

10-2p 

Alxim. 

2.7 

6400 

0.2 



Lucite 

1.18 

2680 







Strength Properties 

Elastic 


Name 

Density 

Tensite Strength Compression 

(N/cm^) (N/cm^) 

Mod-ulus 

Ref. 


Chicken Muscle 
duck tendon 
femur bone 
(domestic fowls) 
Ostrich femur 


59-98 

6370 


[161 

[24] 

6860 [24] 

long direction P 

1.36 MN/cm [241 

(tension) 
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TAELEII. - STRESS -STRAIN COEFFICIENTS FOR 55 PERCENT GRAPHITE 

FIBER-EPOXY MATRIX COMPOSITE 
(All constants to be multiplied hy 10® pel , see Figure T) 
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Figure 22. Water Hanmer Model for Fluid - Solid Impact 
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Figure Moraentu::! Model for Imoact Force due to Fluid Impact 

witli a Flat Plate 


«7 




0.0 SEC 

BB3E3PB3EBBEBEBBBBBBBBBBBBBBBBBBBBBBBBBI 


B 

BS S S S 

BFFFFFFFPS S 

BFFFFFF^FFFFFS 

BFFFFPFPFPFFPFFS 

BFFFFPFFFPFFFFFS 

BFFPFFFFFPFPFFFFFS 

BPPFFFPPFFFFPFFFFS 

BFFFFFFFFFFFFFFFFS 

B FPFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFS 

3PFFFFFPPFFPFFFFFS 

BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFS 

BFFt. f FFPFFPFFFS 

BFFFFFFFFFFFFS 

BFFFFFFFFS S 

BS S S S 

0.30000E-03 SEC 
B3BBBBB3BBBBEBBBEEEEEEE 
BFFFFFFPPFFFFPPPFS S 
8PFPFFFFFFFFFFFS 
BFFFFFFFPFFPFFFFFS 
BFFFPFFFFFFFPFPPPS 
BFFFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFFFFFS 
BFFPf FF?F?FF?''F?FS 
BFFFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFFFFFS 

bffffppfpfffffffps 

BFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFS 
BFFFFFFFFS S 
BS S S S 
3 
3 
3 



U.15000S--03 SEC 

BEBBBBBBBEBBBBBEBBEBBBEBBBBBBEEEEEEEEBE 
BFFFFFFFFS 
BFFFFFPFFFFS S 
BFFFFFFFFFFFFS 
BFFFFFFFFFFFFFFS 


3FFFFPFPPFFFFFPFFS 

BFFFFFFFFFFFFFFFFS 

BFFPFFFFFPPFFFFFPS 

BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFS 


Diansiterj 1.8 cm 

V ; - 10000 cffl/sec 

oy ’ ' 

V ; 9*5^ poise 


BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFS 

BFFFFFFFFFFS S 

BFFFFFFS S 

BS S S 

B 


0.44999S-03 SEC 

BBEBEBEBBBEBPBPEBBEBBBEBBBBBBBBEEEEEEBE 
.BFFFFFFFFFFFFFFFFFFFFS S S S 
BFFFFFFFFFFFFFFFFFFS 

Ibfffpffffffffffffffs -=r.- 

BFFPFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFFFS 

BFFPFFPFFFFFFFFPFFFS ‘ '■ ‘ ' 

BFFFFFFFFFFFFFFFFFFS 
SFFFFPFFFPFFFFFFFPFS 
BFFFFFFFFFFFFFFFFS -- - 

BFFFFFFFFFFFFFFFFS . V;. ' - 

B^FFFFFFFFFFFFFS 

BFFFFFFFFFFFFS ' 

BFFFFFFFFFFS S 
BFFFFFFS S 
BS S S 
B 

B . ■ . 

B ■■ ‘ ' ‘ 

B 


Figure 26a 


RODUCii'U.tt V 


0.59999E-03 SEC 

BBBB5EB3BBBBBEEEBEEEEEEBBBBBBEEEEEBBBBE 
BFPFFFFEFFFPPFFPFPFFFS S S S S S 

SFF^FFFFFFFF’FFFPFPFFS 
3FFEFFFFFPFFFFFFFFFFFS 
BPFFFFFFFFFFPPfFFPFPPS 
BFFFFFFFPFFFFPFFFFFFFS 
BFFFPFFFFFFPFFFFFFFS 
3 pwppw p ffFFFFFFFFFFS 
BFPFFFFFFFFFFFFFFS 
3 F F? ? F F F r ? ?F ^’^FPFS 
BFFFPFFFFFFFFFFFPS 
9FFFFFFPFFPFFFFS 
BFFFFFFFFFFFFS 
3FPFFFFS S S 
3S S S 
B 
3 
S 
B 
3 


0.74998E-03 SEC 

BBBBB5EBBBBEEBEEEBEEEBBBBBBBEEEB EE BEEBE 
BFFFFFFFFFFFFFFFFFFFFFFS S S S S S S 

BFFFEFFFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFFFPFFFS -- . 

EFFFFFFFFFf FFFFFFPFFFS 
BFFFFFFFPFFFFFFFFFFFFS 

BFFFFFFFFFFFFFFFFFFFFS — - - 

3FFFFFFFFFFFFFFFFFFFFS 
BFFFFFFF PFFFPFFPFFFS 
BFFF? FFFFFFFFFFFFS 
BFFFFFFFFFFFFFFS 
BFFFFFFFFFFFFS 

BFFFFS ?FS S 

BS S S 

B 

B 

3 

S 

B 

3 


0.89997E-03 SEC 

BBB3BBBBBBB3BBBBEBBBBEBBBBBBEBBBEEEBBBE 
B FFFFFFPPFFFPFPFFFFFFS FFFFS S S S S 
BFFFFFPFFFFFFPFFFFFS S S 
3?FFP?«?FFPFF?FFFS 
3FPPFFFPPPFFFPFFFPFS 

gpppppwppppppppppppg 

BFFFFFFPPPFPFFFFFFFFFS 

Br FPFPFFFFFFPFFFFPFPFS 

3FPFFPFFPPFPPFFPPPFS 

3FFFFFFFFFFFFFFFFS 

BFFFFFFFFFFFFS S S 

BPFFFFFFFS S 

3S S S S 

B 

3 

3 

3 

a 

3 

3 


0.10400S-02 SEC , 

BEEEB3B3BBB3BBBBBBBBBBBE3BBBBBEBBBBBBBB 
BFFF i’FFFFFFFFFFFFFFFFS S FFFFFFS S S S 
BFFFFFFFFFFFFFPFFPFS S S S 

BFFFFFFFFFFPFFFFFFPS 
3FFFFFFFFFFFFPFFPFFS 

BFFFFPFFFFFf PFPFFFFS , 

BFFFFFFFFFFFFFFFFFFFFS 

3PFFFFFFFFFFFFFFFFFS 

■BFFPFFFFFFFFFFFFFFFS 

BPFFFF? FFFFFFFFFFFFS 

BFFFFFFFFFFFFS S S 

BFFFFFFFFS S . 

BS S S S 

B 

B 

B 

B 

B 

E 

3 


FiffbUrc' ?.6h 





t CO^ WTI yi‘ 
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Impact Veol c i t v= 1000 cm/sec. 
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Figure 31* Time Sequence of 
Velocity Vectors for Normal 
Impact of a Sphere with a Rigid 
Wall (conditions same as in 
Figure 30 ) 
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Figure 32. Velocity Vectors for No.Tnal Impact of a Sphere, (Diameter 
1.8 cm. Velocity 1000 cm/s., Kinematic Viscosity 0.018 
poise) Showing Jevelomenl' of Ediy Flow 
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Figure 33* Radial Velocity Distribution along Wall after 
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Figure 3U. Normal Velocity Distribution along the Lnpact Axis 
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Figure 37 . Pressure Distribution along Vail after Impact /Diameter; 1.8 ccJ 
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